Phase Space Approach to Solving the Time-independent Schrodinger Equation 
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We propose a method for solving the time independent Schrodinger equation based on the von 
Neumann (vN) lattice of phase space Gaussians. By incorporating periodic boundary conditions 
into the vN lattice [F. Dimler et al., New J. Phys. 11, 105052 (2009)] we solve a longstanding 
problem of convergence of the vN method. This opens the door to tailoring quantum calculations 
to the underlying classical phase space structure while retaining the accuracy of the Fourier grid 
basis. The method has the potential to provide enormous numerical savings as the dimensionality 
increases. In the classical limit the method reaches the remarkable efficiency of 1 basis function per 
1 eigenstate. We illustrate the method for a challenging two-dimensional potential where the FGH 
method breaks down. 

PACS numbers: 2.70.Hm, 2.70.Jn, 3.65.Fd 82.20.Wt 



The formal framework for quantum mechanics is an 
infinite dimensional Hilbert space. In any numerical cal- 
culation, however, a wave function is represented in a 
finite dimensional basis set and therefore the choice of 
basis set determines the accuracy. The optimal basis set 
should combine accuracy and flexibility, allowing a small 
number of basis functions to represent the wave functions 
even in the presence of complex boundary conditions and 
geometry. Unfortunately, these two criteria — accuracy 
and efficiency — are usually in conflict, and globally ac- 
curate methods [H-Q lack the flexibility of local methods 
0-0] ■ For example, in the pseudospectral Fourier grid 
method the wave function is represented by its values on 
a finite number of evenly spaced grid points. Due to the 
Nyquist sampling theorem, this allows for an exact repre- 
sentation of the wave provided the wavefunction is band 
limited with finite support |8HlO|. However, the non-local 
form of the basis functions in momentum space leads to 
limited efficiency. On the other hand, in the von Neu- 
mann basis set 11, T3| each basis function is localized on 
a unit cell of size h in phase space. However, despite the 
formal completeness of the vN basis set 13| , attempts to 
utilize this basis in quantum numerical calculations have 
been plagued with numerical errors 0, Ell- 
in this paper we establish a precise mathematical for- 
malism for the vN basis on a truncated phase space. By 
using periodic boundary conditions in the vN basis, as 
introduced in the seminal work by Dimler et al. 15[ , the 



method achieves Fourier accuracy with Gaussian flexi- 
bility. This allows one to tailor the basis in quantum 
eigenvalue calculations to the underlying classical phase 
space structure, with the potential for enormous numer- 
ical savings. The efficiency of the method relative to the 
Fourier grid rises steeply with dimensionality, defeating 
exponential scaling. In the classical limit the method 
reaches the remarkable efficiency of 1 basis function per 
1 eigenstate. 



herent states" of the form: 
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where n and I are integers. Each basis function is a Gaus- 
sian centered at (na, in phase space. The parameter 
a = tt 2 - controls the FWHM of each Gaussian in x and p 
space. Taking Ax = a, Ap = h/a as the spacing between 
neighboring Gaussians in x and p space respectively, we 
note that AxAp = h so we have exactly one basis func- 
tion per unit cell in phase space. As shown in 13[ this 
implies completeness in the Hilbert space. 

The "complete" vN basis, where n and / run over all 
integers, spans the infinite Hilbert space. In any numeri- 
cal calculation, however, n and I take on a finite number 
of values, producing N Gaussian basis functions {<7;(x)}, 
i = 1...N. Since the size of one vN unit cell is h, the area 
of the truncated vN lattice is given by S* vN = Nh. 

The pseudospectral Fourier method (also known as the 
sine Discrete Variable Representati on [16lp is a widely 
used tool in quantum simulations |l7H20j . In this method 
a function tp(x) that is periodic in L and band limited 
in K = j- can be written in the following form: tfj(x) — 
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jr. The basis functions {9 n (x)} are given by |2l| : 
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which can be shown to be sine functions that are periodic 
on the domain [0,L] @. The set {^(x)} i = 1,..,N 
spans a rectangular shape in phase space with area of 
5 FGH = 2LP = 2L^ = Nh. Thus N unit cells in the 
vN lattice and N grid points in the Fourier method cover 
the same rectangle with an area in phase space of: 



The von Neumann basis set 121 is a subset of the "co 
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(Fig. [TJ. This suggests that N vN basis functions con- 
fined to this area will be equivalent to the Fourier basis 
set. Unfortunately, the attempt to use N Gaussians as 
a basis set for the area in eq.© (Fig. [T]) is unsuccessful, 
a consequence of the Gaussians on the edges protruding 
from the truncated space. However, by combining the 
Gaussian and the Fourier basis functions we can generate 
a "Gaussian-like" basis set that is confined to the trun- 
cated space. We use the basis sets {gi(x)} and {8i(x)} 
to construct a new basis set, {(ji(x)}: 
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FIG. 1: N = 9 coordinate grid points and N = 9 vN unit cells 
(4) span the same area in phase space, 5 = Nh. The vN basis functions 
are Gaussians located at the center of each unit cell. 



for m = 1, N. The new basis set is in some sense, the 
Gaussian functions with periodic boundary conditions. 
We can write eq.(|4]) in matrix notation as: G = QG 
where Gij — gj{xi) By taking the width parameter a — 
we can guarantee that the pvN functions have no 
linear dependence and that the matrix G is invertiblc, 
that is GG^ 1 — 9. The invertibility of G implies that 
both bases span the same space. 

The representation of in the pvN basis set is given 
by: 



lar, if \if>) has an irregular phase space shape we may 
expect that some of the pvN basis functions will ful- 
fill the relation: (jjj\ip) = 0, j = 1,...,M. Due to the 
non-orthogality of the basis we cannot simply eliminate 
the states §j, since the coefficients of gj may include 
contributions from remote basis functions, but we can 
take advantage of the vanishing overlaps by defining a 
bi-orthogonal von Neumann basis (bvN) {bi(x)}. 
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To find the coefficients a m we first define the overlap 
matrix, S: 
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Using the completeness relationship for non-orthogonal 
bases, \ip) can be expressed as 
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Comparing with 
En=i( 5 '~ 1 )'™™(5n|V') and (</, 



cq.© we find that 

Although the periodic von Neumann (pvN) and the 
Fourier methods span the same rectangle in phase space, 
the localized nature of the basis functions in the pvN 
method can lead to significant advantages. In particu- 
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or in matrix notation: B = GS . 
eqEl can be written as 



Inserting eq[9] into 
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By assumption, M of the coefficients are zero, hence in 
order to represent in the bvN basis set we need only 
N' = N~M basis functions. Note that the bvN and pvN 
are bi-orthogonal bases, meaning that each set taken by 
itself is non-orthogonal but they are orthogonal to each 
other. This can be shown easily by: 
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For many practical applications the full knowledge of the 
basis wavefunctions is unnecessary: wc need only the 
value of the basis functions at the sampling points. For 
example the evaluation of Hamiltonian matrix elements 
can be performed explicitly by: 
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and similarly: 
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where _ff FGH = t/FGH + T FGH and the potential and the 
kinetic matrix are given by: Vu « V(xi)8ij and 
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|23j. The eigenvalue problem in a non-orthogonal basis 
set becomes HU = slIE; in the pvN basis set s is given 
by eq. and in the bvN basis set s is given by: 



B*B = S-i&GS- 1 = S- 1 . 



(15) 
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Diagonalization should give accurate results for all wave- 
functions localized to the classically allowed region of the 
rectangle. Note that in the multidimensional implemen- 
tation, the S^ 1 matrix required in Eq.(9) can be con- 
structed separately for each dimension. As a result, the 
computational effort to construct the bvN basis set is 
negligible compared with diagonalizing the Hamiltonian. 

As a numerical test of the pvN basis set we stud- 
ied the standard example of the harmonic oscillator 
V(x) = mu 2 x in units such that m = % = uj = 1. 
We calculated the seventh excited energy using differ- 
ent number of pvN and conventional Gaussian basis _14,3 
functions. In the Gaussian basis set the Hamiltonian 
and the overlap matrices were calculated analytically as: 
H l3 = ( gi \H\ gj ) =J™ oo g* i {x)[-£s + V{x)] gj {x)dx and 
Sij = (9i\9j) = I X 00 9i(x)9](x)dx. The results, shown 
in Fig. UKa), show the superiority of the pvN basis set 
compared to the standard Gaussian basis set. In fact, 
the results obtained with the pvN basis set are exactly 
as accurate as in the Fourier grid method. The kinetic 
energy spectra in Fig. G{b) show that the pvN has a 
perfect quadratic dependence while the vN spectrum is 
highly flawed. 
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FIG. 2: (a) Error in the 7 th eigenvalue of the harmonic oscilla- 
tor as a function of basis set size for vN(dashed) and pvN(solid). 
(b) Kinetic energy spectra using 16 basis functions. vN(dashed), 
pvN (solid). 



both the FGH and bvN basis sets. The Morse parameters 
were taken to be D = 12, m — 6, (3 — 0.5 and H = 1. For 
FGH, 100 grid points between [—1.6,20.1] were required 
to get 4 digits of accuracy in energy for all 24 bound 
states. By using the bvN basis functions (constructed 
from 10x10 vN functions with a = 0.5) we obtain the 
same 4 digit accuracy with only 48 basis functions. This 
is demonstrated graphically in Fig. (a). The figure 
shows the phase space representation of 100 evenly grid 
points. Although it requires 100 pvN basis functions to 
span this area in phase space, due to the flexibility of the 
bvN basis set we can suffice with just the basis functions 
in the classically allowed region (magenta squares). 
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FIG. 3: (a). Phase space area spanned in the bvN method (ma- 
genta) and in the pvN (or FGH) method (full rectangle) for Morse, 
(b) Efficiency ratio (defined as number of basis functions per con- 
verged eigenstates) of the pvN (solid) and FGH (dashed) methods 
for the Morse potential as function of h. 



In the bvN basis set we are able to remove some of the 
basis functions and construct lower dimensional H hvN 
and 5* bvN matrices without losing accuracy. In order to 
test this claim, we calculated numerically the eigenener- 
gies of the Morse oscillator V(x) = D(l — e~ l3x ) 2 by using 



The ability to localize a bvN function at a specific point 
in phase space results in the remarkable concept of 1 ba- 
sis function per 1 eigenstate. This means that in order 
to calculate N eigenenergies we need only N basis func- 
tions. Obviously, such one per one efficiency, if reachable, 
will be the ideal efficiency for any basis set. In order to 
test the ability of the bvN method to reach the ideal ef- 
ficiency we examined the Morse potential and looked for 




FIG. 4: The triangle potential: a two dimensional test case for the 
bvN method. 



the smallest basis that provides 12 digits of accuracy for 
all the eigenvalues up to E = 11.25. The bvN method 
indeed tends to the ideal efficiency in the classical limit 
H — > (Fig. [5b). This remarkable result is unique for 
methods based on phase space localization 24|. 

The true power of the method is in the application 
to higher dimensional systems. As an illustration, con- 



sider the potential V(r, 6) = (1 — cxp(— a(6)r 2 )) 2 where 
a = ((1 - cos(3<9))/4) 2 +0.05. This 3-fold symmetric po- 
tential (Fig. (01), which is a realistic model for a system 
of three identical particles and fixed hyperradius, is quite 
challenging for the FGH method. Taking m = 96, H = 1 
gives 760 states below E = 0.996. In order to get two 
digits of accuracy for all those states one needs ~ 11000 
FGH grid points, while with the bvN basis set conver- 
gence is achieved with only 1500 basis functions. For 
higher accuracy (four digits), the FGH breaks down com- 
pletely while the bvN method requires fewer than 3000 
basis functions (Fig.fl5h.-b)). Figure fl5b)) shows again 
that as h — > the efficiency tends to 1 basis function 
per 1 eigenstate (because of the size of the calculations 
we consider just 3 digits of accuracy). In contrast, the 
FGH efficiency as H —> is determined by the ratio be- 
tween the classical phase space and the box that contains 
it, which we calculate to be ~ 10 for this system using 
Monte Carlo integration. Note that in going from 1-d to 
2-d the savings provided by the bvN relative to the FGH 
method has gone from 2 to 7-10 for qualitatively simi- 
larly potentials. This suggests that the relative efficiency 
of the bvN method increases rapidly with dimension. 

To explore the scaling with dimensionality more fully, 
consider a harmonic oscillator with 1-d classical phase 
space volume v up to energy E. For the D-dimensional 
oscillator, the total phase space volume up to energy 
E is V = v /D\ In the classical limit, the total num- 
ber of states is determined by V/h D and therefore in 
this limit the efficiency of pvN relative to FGH is deter- 
mined by the ratio of the phase space volumes spanned. 
Defining a to be the area of the box surrounding the 
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FIG. 5: Triangle potential results for bvN(solid) and FGH(dashed) 
(a) The calculated highest eigenenergy as a function of basis set 
size N. (b) The accuracy of the calculated highest eigenenergy as a 
function of basis set size N . (c) Efficiency ratio of the pvN (A) and 
FGH (O) methods as a function of h. The A (pvN) and • (FGH) 
signify that the value is an approximation to the h — > value, given 
by the ratio between the size of the phase space spanned by the 
basis and the classical phase space. 



1-d oscillator phase space, the volume of the box sur- 
rounding the D-dimensional phase space is A = a D and 
the ratio of phase space volumes is S = V/A = s D /D\ 
where s — v/a — 7r/4 for the harmonic oscillator. For 
the Morse, Coulomb and other chemically relevant po- 
tentials, the 1-d ratio s < ir/4 and the D-dimensional 
phase space volume scales more slowly than v D /D\ [22l |; 
these effects combine so that the relative efficiency of the 
pvN method rises steeply with dimension. As a result of 
the D\ in the expression for V, the method remarkably 
defeats exponential scaling. A more detailed analysis [22j 
shows that for D 3> v/h — g the method scales polyno- 
mially: V = D 9 /g\ 

Work in progress includes application to vibrational 
eigenvalue calculations for realistic polyatomic molecules, 
electronic eigenvalues for multiclcctron atoms and exten- 
sion of the approach to the time-dependent Schrodinger 
equation. 

This work was supported by the Israel Science Founda- 
tion and made possible in part by the historic generosity 
of the Harold Pcrlman family. We thank Bill Poirier for 
helpful discussions. 



[1] R.Kosloff in Numerical Grid Methods and their Applica- 
tion to Schrddinger's Equation ed. C. Cerjan (Kluwer, 
Boston, 1993). 

[2] C. C. Marston and G. G. Balint-Kurti, J.Chem.Phys. 6, 

3571 (1989). 
[3] G. W. Wei, J.Phys.B. 33, 343 (2000). 
[4] M.J.Davis and E J.Heller, J.Chem. Phys. 71, 3383 

(1979). 

[5] I. P. Hamilton and J. C. Light, J.Chem. Phys. 84, 306 
(1986). 

[6] Z.Bacic, R.M Whitnell, D.Brown and J.C.Light, Comp. 

Phys. Comm. 51, 35 (1988). 
[7] S. Garashchuk and J. C. Light, J.Chem. Phys. 114, 3929 

(2001). 



5 



[8] E. T. Whittaker, Proc. R. Soc. Edinburgh 35, 181 (1915). 
[9] H. Nyquist, Trans. AIEE 1 47, 617 (1928). 
[10] C. E. Shannon, Proc. IRE 37, 10 (1949). 
[11] S. Fechner, F. Dimler, T. Brixner, G. Gerber and D. 

J.Tannor, Opt. Express 15, 15389 (2007). 
[12] J.von Neumann, Math. Ann. 104, 570 (1931). 
[13] A.M Perelomov, Theor. Math.Phys 11, 156 (1971). 
[14] B. Poirier and A.Salam, J.Chem. Phys. 121, 1690 (2004). 
[15] F. Dimler, S. Fechner, A. Rodenberg, T. Brixner, and D. 

J.Tannor, New J. of Phys. 11, 105052 (2009). 
[16] D. T. Colbert and William H. Miller, J.Chem. Phys. 96, 
1982 (1992). 

[17] J. Dai and J. C. Light, J.Chem. Phys. 107, 1676 (1997) 
[18] A. J. H. M. Meijer and E. M. Goldfield, J.Chem. Phys. 

108, 5404 (1998) 
[19] X. T. Wu, A. B. McCoy and E. F. Hayes J.Chem. Phys. 

110, 2354 (1999) (2002) 
[20] J. H. Baraban, A.R. Beck, A.H. Steeves, J.F. Stanton 

and R. W. Field, J.Chem. Phys. 134,244311 (2011) 
[21] D. J. Tannor, Introduction to Quantum Mechanics: A 

Time-dependent Perspective (University Science Books, 

Sausalito, 2007), eq.11.163. 
[22] See attached supplementary material. 
[23] Ref.[lU eq.11.172. 

[24] R. Lombardini and B. Poirier, Phys Rev E. 74, 036705 
(2006) 



SUPPLEMENTARY MATERIAL 



(N 

o 



43 

a: 



> 

(N 
(N 

o 

(N 



X 



I. THE DIRICHLET OR PERIODIC SINC 
FUNCTION 

In Eqn. (2) of the paper we introduced the pseudospec- 
tral functions underlying the FGH/sinc DVR method as 
a sum of exponentials: 



»(*) 




i2irj 

exp ( ——(x-x n ) . (1) 



In this section we will show that these pseudospectral 
functions 9{x) are periodic sine functions, proportional 
to the Dirichlct functions. We will first replace the index 
j with k = j + N/2 — 1 and define A = 27T ( x ~ Xn ') (f or 
simplicity we suppress the subscript n on A) to obtain: 
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■. exp (iA(k- N/2 + 1)) 



N-l 

lK^ +1 )) g (exp(iA)A (2) 



fc=0 



For yl = 27rm, where m is an integer, 6 n (x) = \ jf] for 



A 7^ 27rm we use the expression for a geometric series 



fe=0 



to get: 



exp (i^(-f + 1)) 1 - exp(iAN) 
\J~LN 1 - exp(iA) 

exp(^) sin(fj 



V^LiV sin(f) 



where -Dat(A) is the Dirichlet function. 



II. SCALING OF THE PVN METHOD WITH 
DIMENSIONALITY 



ciency of the pvN relative to the FGH method is (the 
inverse of) the ratio of the classical phase space volume 
to that of the enclosing box. For harmonic potentials, 
this relative efficiency increases faster than exponentially 
with dimension D. For anharmonic potentials, two fac- 
tors contribute to make the relative efficiency even higher 
than in the harmonic case. 

In Section [II CI we go beyond the semiclassical formu- 
lae. Starting with a fully quantum formula for the num- 
ber of states of a D-dimensional isotropic harmonic oscil- 
lator, we show that the scaling with D is sub-exponential, 
leading in the limit of large D to polynomial scaling, 
G = D'/gl 



Phase space volume V of the D-dimensional 
oscillator 



The phase space of the D-dimensional isotropic har- 
monic oscillator up to energy E, if + ^fi ( w « — 
rnj = 1) is a sphere with radius R = \J2E and volume 

V = —/j< — = — D i ■ I his gives volume v = nR 
1-dimension and volume 



D\ 



(3) 



in D-dimensions (Fig(JTJ). 

In 1-d, the relation between states and volume is g = \ . 
Thus, we obtain a formula for the scaling of the total 
number of states similar to that for the phase space vol- 
ume: 



G 



V_ 

hP 



D 



fj_ 

D\ 



(4) 



Because of the D\ in Eq.@, G increases slower than ex- 
ponentially. In fact according to Eq.Q G decreases for 
D > g; however as we will show in Section HTCl this effect 
is an artifact of using the semiclassical formula outside 
of its domain of applicability. 

80 1 



In this section we discuss the scaling of the pvN method 
with the dimensionality of the system. We begin in Sec- 
tion III Al by discussing the scaling of the classical phase 
space volume V up to energy E of a set of D isotropic 
harmonic oscillators. We show that V increases more 
slowly than exponential with D. 

Since the number of states G is given semiclassically by 
G = -Xy, and since in the classical limit the pvN method 
requires only 1 basis function per eigenstate, we find that 
the scaling of the pvN method is slower than exponential 
with D. 

In Section III Bt we note that semiclassically the effi- 
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FIG. 1: Log of the volume of the classical phase space below E = 8 
as a function of dimensionality (solid) compared with exponential 
scaling(dashed) . 
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B. Ratio between phase space volume and box 

1. Harmonic potentials 

The efficiency of the pvN method relative to the FGH 
in the classical limit is (the inverse of) the ratio of the 
classical phase space area and the enclosing box. For 
a 1-d harmonic oscillator, this is the ratio between the 
area of a circle and the square surrounding it, given 
ky : s = (Ufp = <P ^ n -D-dimensions, this ratio is 

S = jyr = ■ Because of the D\, the relative ef- 

ficiency (the inverse of S) increases faster than exponen- 
tially with dimension D. 
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FIG. 3: Log of the phase space volume V vs. D for a D- 
dimensional sum of Morse potentials, calculated using Monte Carlo 
integration (dots) and using Eq.Q (solid). Exponential scaling is 
demonstrated by the dashed line. 



2. 1-d anharmonic potentials 



C. Number of states for D-dimensional potential 



For the harmonic oscillator, the packing parameter 
s = j w 0.7854 is relatively close to 1. However most 
chemically relevant potentials have classical phase space 
shapes that are packed in a much less compact way 
(Fig.®). For example, for the Morse oscillator shown 
in Fig[2][b) s — 0.4638; for the Coulomb potential shown 
in FigHJc) s = 0.0321. 




FIG. 2: Phase space area of (a) harmonic, (b) Morse and (c) 
Coulomb potentials. 



3. D-d anharmonic potentials 

For anharmonic potentials, not only is s < ?, but 
the multidimensional classical phase space volume scales 
more slowly than in Eq. ([3|) . This can be understood by 
visualizing figures such as Figs. EJb)-(c) in multidimen- 
sions, and is confirmed numerically using Monte Carlo 
integration (see Fig.©); in Section Til CI we rationalize 
this behavior in a different way using state counting ar- 
guments. As a consequence, the ratio of phase space 

volumes Saho is even smaller than iSho = jyr (and the 
efficiency of pvN is therefore larger): 



A HO 

Dl 
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Dl 



> Saho- 



(5) 



Since the harmonic efficiency (the inverse of Sho) in- 
creases faster than exponentially with dimension D, 
clearly the anharmonic efficiency (the inverse of Saho) 
does as well. 



In Section III Al we showed that for the D-dimcnsional 
isotropic harmonic oscillator, in the semiclassical limit, 
the number of states below energy E is given by G — 
g D /D\ We now reexamine this result starting from a fully 
quantum treatment. For a 1-d oscillator the number of 
states below energy E is g « For a D-dimensional 
isotropic harmonic oscillator the number of states below 
E is not g D , since some combinations give energy above 
E. Computing G is identical to the number of ways to 
distribute no more than g indistinguishable balls (quanta 
of energy) into D distinguishable boxes (oscillators) . The 
answer is given by [l|: 



G = J2 



i + D-l (g + D)\ 



i\{D-l)\ 



g\D\ 



(G) 



We now draw attention to two points about Eqn. [6] 

1) Consider the case D = 2 with g 2. Then Eqn. [6] 

becomes 



G 



(g + 2)(g + l) _ g 2 



(7) 



The factor g 2 in the numerator has its origin in the ten- 
sor product Hilbert space of the two truncated harmonic 
oscillators. The factor of 2 in the denominator comes 
from the fact that about half of the states in the tensor 
product Hilbert space have energy E\ + E^ > E. For 
a D-dimensional system the fraction of tensor product 
states with Yjf =l E.i < E is approximately 1/iV! 

In the case of an anharmonic potential such as the 
Morse oscillator, the gaps between energy levels are not 
equal: the higher energy states are spaced more closely 
together. As a consequence, for a larger fraction of com- 
binations than in the harmonic case HfL 1 Ei > E and the 
scaling of G is slower than in Eq.(g]) (Fig. ([3])). 
2) If we assume g 3> D, Eq.© becomes 



D 



G 



fj_ 
Dl 



(8) 



3 



which is the same result as Eq.Q. However, in the op- 
posite limit, D 3> g, we obtain 



G 



D3 



(9) 



showing that for large enough D the scaling D is actu- 
ally polynomial. Equations [SHU are plotted in Fig. 2] 
The true scaling with D is clearly subexponential. Equa- 
tion [5] is seen to be a valid approximation for g D, 
although the turnover at large D is an artifact of using 
the semiclassical formula outside of its range of validity. 
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FIG. 4: Log of the number of states G that have enegy below E = 
ghuj, (g = 30) vs. D for a D-dimensional harmonic potential. The 
figure shows the exact expression (red) , the approximations for g S> 
D (blue), for D ^> g (green) and a comparison with exponential 
scaling (dashed). 
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